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Abstract 

The  state  estimation  problem  for  bilinear 
stochastic  systems  evolving  on  spheres  Is  con- 
sidered. The  problem  Is  motivated  by  some  appli- 
cations involving  rotational  processes  in  three 
dimensions . Then  the  theory  of  harmonic  analysis 
on  spheres  is  used  to  define  assumed  density 
approximations  which  result  in  Implementable 
suboptlmal  estimators  for  the  state  of  the  bilinear 
system. 

1 . Introduction 

Fourier  series  analysis  has  bean  applied 
In  several  recent  studies  [l]-[4]  to  estimation 
problems  for  stochastic  processes  evolving  on 

the  circle  . Wlllsky  [4]  used  Fourier  series 
methods  to  define  "assumed  density"  approxima- 
tions for  certain  phase  tracking  and  demodulation 
problems.  In  feet,  a system  designed  using  these 
techniques  performed  better  than  other  estimators , 
including  an  optimal  phase-lock  loop. 

In  this  paper  we  study  bilinear  systems 
evolving  on  spIMraa  (the  more  general  case  of 
systems  evolving  on  compact  Lie  groups  or  homo- 
geneous spaces  la  discussed  In  [13]).  The  optimal 
estimator  Is  In  general  Infinite  dimensional  [7], 
and  our  approach  to  the  design  of  suboptlmal 
estimators  Is  a generalization  of  that  of  Wlllsky 
[4] , whose  work  is  described  in  Section  3 . The 
basic  approach  Involves  the  definition  of  an 
"assumed  density"  form  for  the  conditional  density 
of  the  system  stata  at  time  t given  observations 
up  to  time  t.  These  densities  are  defined  via  the 
techniques  of  harmonic  analysis  on  spheres  [S] , 
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[6l  (which  generalize  the  Fourier  series  on  the  Lie 

group  sS«  Our  method  differs  from  most  previous 
assumed  density  approximations  in  that  our  ap- 
proximation Is  defined  on  the  approprlata  compact 
manifold  (as  opposed  to  Gaussian  approximations . 

for  example,  which  are  defined  on  R*'  [7l). 

In  Section  2 we  discuss  some  general 
properties  of  stochastic  bilinear  systems  and  pre- 
sent an  example  Involving  satellite  tracking. 

The  use  of  harmonic  analysis  will  be  moUvalad  by 
the  phase  tracking  example  of  Wlllsky  [4]  In  Sec- 
tion 3.  In  Section  4,  wc  discuss  the  optimal 

2 

estimation  problem  on  the  2-sphere  S and  the 
application  of  harmonic  analysis  to  the  design  of 
suboptlmal  estimators . 

2 . Stochastic  Bilinear  Systems 

* The  basic  stochastic  bilinear  systam  con- 
sidered here  Is  described  by  the  Ito  stochastic 
differential  equation  [4],  [9],  [lOl,  [12],  [13], 
[14l-[16l,  [17].  [18]. 


dx(t)  « {[Ajj  + Qy(t)AjA^]dt 

+ (0]x(t)  (1) 

i-l 

where  x Is  an  n-vector,  A^  are  n x n matrices. 

Is  the  (l,J)th  element  of  Q.  and  w Is  a Browni- 
an motion  (Wiener)  process  with  strength  Q(t) 
such  that  E[w(t)w’(s)]  ■ J^"*^"^*'*^Q(T>dT  . Fol- 
lowing the  notation  of  [8]-[10],  [11],  we  define 
X “ (Ag , Aj , . . . , A^  to  be  the  Lie  algebra  con- 
taining these  matrices . The  corresponding  .000*^ 
nected  matrix  Lie  group  Is  denoted'by  0-{ext>X]Q. 
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Thea  the  a^utlon  x(t)  of  (1)  will  evolve  on  the 
homogeneous  space  G*x(0)  (l.e.,  x(t)CGx(0)  for 
all  txO)  in  the  mean-square  sense  and  almost 
surely  [isl.  lid. 

Associated  with  the  Ito  equation  (1)  is  a 
sequence  of  equations  for  the  moments  of  the  state 
x(t),  first  derived  by  Brockett  (9  1,  [10],  [14],  if 

N(n,p)  denotes  the  binomial  coefficient  ^ , 

then  given  an  n-vector  x,  we  define  x^**^  to  be  the 
N(n,p)~ vector  with  componenU 


(ordered  lexicographically) . These  components 
are  the  monomials  ( or  homogeneous  polynomials) 
of  degree  p In  Xj x^,  the  components  of  x, 

scaled  so  that  I|x|^  - Bx^**^!!.  Given  an  mxn 

mafrlx  A.  wo  denote  by  A^^^  the  unique  matrix 
which  veriflos 


y- Ax  ^y' 


,Ip1  . j^lpl  j,[pl  , 


It  U clear  that  11  x saUsfies  ttie  linear  dlfferenUal 
equation 


x(t)  - Ax(t) 


(4) 


The  observation  model  consldored  In  this 
paper  conslsU  of  linear  observations  of  the  state 
corrupted  by  additive  urhlte  noise,  or 

ds(t)  - H(t)x(t)dt  ♦ R^(l)dv(t)  (7) 

where  s Is  a p-vector,  R>0,  and  v is  a Wiener 
process.  This  bilinear  system-linear  observation 
model  is  useful  In  the  study  of  certain  practical 
problems,  as  lllusbated  by  the  following  example. 

Consider  a satellite  In  circular  orbit  about' 
some  celesttal  body.  Because  of  a variety  of 
effects  Including  anomalies  In  the  grevltetlonal 
field  of  the  body,  effects  of  the  gravitational 
fields  of  nearby  bodies,  and  the  effects  of  solar 
pressure,  the  orbit  of  the  satellite  Is  perturbed, 
in  this  case,  the  position  x of  the  satellite  can 
be  described  by  the  simplified  bilinear  model  [18l 


(9) 

form  a basis  for  so(3) . the  Us  algebra  of  3x3  skew- 
symmetric  matrices . Also.  ^ are  the  components 


of  the  nominal  angular  velocity  and  are  the 

components  of  a Wiener  process  with  strength  Q(t). 
If  ^'(0)x(0)l*l.  then  E(x‘(t)x(t)l-1  for  ell  t;  thus 

X evolves  on  the  2-sphere  (the  sanw  statement 
can  be  made  almost  surely  [is], [16]).  We  note 
that  the  assumption  In  (8)  that  ttM  perturbations 
In  the  angular  velocity  are  white  Is  a slmpUflca- 
Uon  [131.  but  the  simplified  model  (8)  can  lead  to 
simple  but  accurate  on-line  tracking  scbssMS.  If 
we  are  then  given  the  noisy  observations  (7)  of 
the  satellite  position,  the  problem  Is  to  estbsate 

x(t)  given  the  past  observations  a*  ^ (s(s)  .Osset) 
3.  A Phase  Tracking  Problem  on  8^ 


then  x*^^'  also  saUsfles  a linear  differential  equa 
tlon 


We  regard  this  as  the  definition  of  A^j . which  Is 

the  Infinitesimal  version  of  A^^^.  In  fact,  A^j  can 
be  easily  computed  from  A [17].  Brockett  has 

shown  that  If  x saUsfles  (1),  Oian  x^^^  saUsfies 
the  Ito  equaUon 


We  first  discuss  a phase  tracking  problem 
studied  by  Bucy,  et  al,  [ll  and  Wlllsky  [4],  In 
which  the  phase  8 and  the  observaUon  s are 
described  by 


de(t)  ■ a dt  q*(t)dw(fi.  t(0)>B 


As  we  shall  see  In  the  sequel,  the  sequence  of 
moment  equaUons  Is  a valuable  tool  In  the  study  of 
state  estimation. 


i 


i 


j 

■1 

i 

■ ■i 
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whore  v and  w an  Independent  standard  Brownian 
motion  processes  Independent  of  the  random  Initial 
phase  t..  We  wish  to  estimate  PCt)  mod2iT  given 

s . and  we  take  as  our  optimal  estimation  criterion 
the  minimization  of 

E[(l-cos(e(t)  - V(t))|z(s),  Oasst]  (12) 

where  this  notation  denotes  the  conditional  expec~ 
tatlon  given  the  e-fleld  o[x(s).  Oasst)  generated 
by  the  observation  process  up  to  time  t. 

Noting  that  we  are  essentially  tracking  a 

point  on  the  unit  circle  S^,  we  reformulate  the 
problem  In  Cartesian  coordinates . Let  - sin  l(t) , 

X2  “ cos  S(t) . Then 

Xj(^pq(t)dt/2  •gdtfq*(t)dw(^jxj(^ 

Xj(tj  [-(•^.dt+q*(t)dw(t))  -q(t)dt/2  JpCjWj 

03) 

dz(t)  ■ Xj(t)dt  + r^(0dv(t)  (14) 

which  are  of  the  bilinear  process  - linear  measure- 
ment type  discussed  In  Section  2 , 


In  Cartesian  coordinates  our  estimation 
problem  Is  to  choose  an  estimate  (x^(t),]^(t))  on 

the  unit  circle.  If  we  use  0>e  constrained  least 
squares  criterion 

I-»E[(Xj(t)-Xj(t))*+(Xj(0->4(i))*  |z(s).0asstl  (IS) 

subject  to  Xj*(t)+Sc!j*(t)-l,  or  equivalently  subject 

to  Xj(t)  » sin  ?^t).  x^(t)  • cos  IKt)  our  criterion 

reduces  to  the  criterion  (12)  used  In  [1]  and  [4). 
One  can  show  [1],  [4l  that 


or 

Kt)  - tan'‘^(Xj(t|t)/x2(t|t))  (17) 

where  the  conditional  expaetttlon  Is  denoted  by 
the  equivalent  notations 


x^(t|t)AB(x^(t)|z(s),  0ss<0^EV(t)l  . (18) 

As  discussed  In  [1]  and  [4] , the  optimal 
(constrained  least-squares)  filter  Is  described  as 
follows . The  conditional  probability  density  of 
• given  fz(s),  Osset)  can  be  expanded  In  the 
Fourier  series 

P(*.t)  - £ c^W**"* 


c^(t)m^Ele'“**^'^  |z(s).0ssst]  Ab  (t)-ta  ft).  (20) 

n 4 rr  n n 


Then  the  optimal  filter  Is  given  by 

2 

dc  (t)  ■ -line  ♦ ^ q(t)lc  (t}dt 
n c r . n 

c ,{t)-c  ,(t)  ^ 

+ ■ + 2«c^(t)Im(Cj(t))  (21) 

T(t)  » tan”*(a^(t)/bj(i))  , (22) 

vrhere  dv(Q  > dz(t)-«-2rIffl(c^(4)dt  . 

The  structure  of  die  optimal  filter  Is  that  of 
an  Infinite  bank  of  filters,  but  the  equation  for  c 

n 

Is  coupled  only  to  the  filters  for  c, . c , , and 

i ll*A 

c ^ . This  fact  plays  an  Important  part  In  the 

lir  1 

following  approximation.  In  order  to  construct  a 
finlte-dlmenslonal  suboptlmal  filter,  we  vrlsh  to 
approximate  the  conditional  deirslty  (19)  by  a den- 
sity determined  by  a finite  set  of  parameters. 
Following  Wlllsky  [4l.  we  assume  that  p(9,Q  Is  a 
folded  normal  density  v/lth  mode  IKt)  and  "varlBitce'’ 

Y(t): 

- F(i;  T(0.  v(t))  . (23) 


The  folded  normal  density  Is  the  solutlan  of  the 
standard  diffusion  equation  on  the  circle  (l.e. , it 

is  the  density  for  Brownian  motion  processes) 

and  is  as  Important  a density  on  as  the  normal 

Is  on  this  point  will  be  dlscuasad  in  more  de- 
tail In  the  next  section.  In  this  ease,  if  c^  lias 

been  computed  and  If  p(*,t)  satisfies  (23)  then 
c^^  can  be  computed  (for  any  N)  from  the  equation 


c 


N+1 


. ,N(N+1)^  (N+1) 


(24) 


Thus  the  bank  of  fitters  described  by  (21)  can  be 
truncated  by  approximating  ^ f>F  (24)  and  sub- 
stituting fills  approximation  into  the  equation  for 
C|^ . This  was  done  for  N>1  in  [4] , and  the  result- 
ing Fourier  coefficient  filter  (FCF)  was  compared 
to  a phase-lock  loop  end  to  the  Gustafson-Speyer 
"state-dependent  noise  filter"  (SDNF)  (l9l.  The 
FCF  performed  consistently  better  than  the  other 
systems,  although  the  SDNF  performance  was  quite 
close. 
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4.  Esttmation  on  S*  ■ P®*"*  ®**  *®  *®™“  ®*  ••»•  Cartoalan 

coordinates  x » (Xj,X2.Xj)  or  the  polar  coordinates 

This  section  Is  devoted  to  the  study  of  the  {t.tfi  (see  the  Appendix,  in  which  hacinonle  analy- 
estlmatlon  problem  tor  the  system  with  the  bilinear  . .2  . . 

state  equation  (1)  with  the  linear  measurement  (7),  ®"  ® *•  summarised).  We  first  review  the 

where  the  matrices  {A,;  1-1. 2. 3)  are  given  by  the  notions  of  Brownian  motion  and  Gaussian  densities 
i on  homogeneous  spaces . «^ch  have  received  much 

skew-symmetric  matrices  IR^;  1-1. 2, 3)  of  (9).  a ttenUon  in  the  literature  (see  [S]. (21]).  Yosida 

This  system  is  moUvatcd  by  the  satellite  tracking  ^ ®*  ® Bf®'»«lan 

problem  of  SecUon  2;  as  discussed  there,  the  state  P"®®*”  °"  »^R‘«»‘«nnlan  homogeneous 

x(t)  evolves  on  the  2-sphere  S*.  ‘P«®®  ” (such  as  S*)  with  respect  to  the  Rleman-  ' 

nlan  measure  is  the  fundamental  solution  of 


The  estimation  criterion  which  will  be  used 
for  this  problem  is  the  constrained  least-squares 
criterion,  which  is  analogous  to  the  criterion  dis- 
cussed in  Section  3 for  the  phase  estimation  prob- 
lem. That  is.  we  seek  R^(t|t)  which  minimizes 

I-  E((x(t)  - ?(t|t))'  (x(t)  - x(t|t))  |sjl  (25) 

* « 
subject  to  Qm  S constraint  Hx(t|t)||  - 

x(t|t)'x(t|t)  - 1.  It  is  easily  shown  (l3l  that  the 
optimal  estimate  la 

x(t|t)-x(t|t)/nx(t|t)||  (26) 


»p(x.  t)/3t  - G*  p(x.  t)  - 0 (29) 

where  G*  is  the  formal  adjoint  of  a differential 
operator  expressible  In  local  coordinates  as 

with  constant  f and  Q - Q*  a 0.  In  partleular.  if  G 

is  the  Laplace-Betmml  operator  (Laplaclan)  Sg2  on 

(see  (A.  4))  which  is  self-adjoint  [51.  the  funda- 
mental solution  of 


erhore  denotes  conditional  expectation  given  the 
e-fleld  generated  by  the  past  observations  (see 
(18)).  Thus  wo  must  compute  the  conditional  ex- 
pectation of  the  state  x(t)  given  the  past  observa- 
tions z^  > (s(s)  .Osset]. 


»p(9.e.fi/»t  - YAp(t.g.l)  - 0 (31) 

where  T > 0.  is  a Brownian  motion  an  S^.  Aeoocd- 
Ing  to  tS].  the  fundamental  solution  of  (31)  is 
given  by 


> 


The  equations  for  computing  the  conditional 
expectation  can  be  derived  from  the  general  non- 
linear filtering  equation  (7 . p.  184]  and  the  moment 
equation  (6) . The  resultant  equations  are 


♦ (E^x^P^(t)x'(t)]-E*(x^V)]E*(x*(t)]  ]H’(t)R"^t)dv(t) 

(27) 

dv(t)  - ds(t)  - H(t)x(t|t)dt  . (28) 


The  structure  of  these  equations  is  quits 
similar  to  that  of  (21)— l.e. . the  estimator  consists 
of  an  infinite  bank  of  filters,  and  the  filter  for  the 

p*^  moment  is  oouplad  only  to  fiiose  for  the  first 

and  (p*^l)*^  moments.  Therefore,  we  are  led  to 
the  design  of  suboptimal  estimators.  Motivated  by 
the  success  of  Bucy  and  Willsky's  phase  tracking 

example  evolving  on  sV  we  would  like  to  design 

2 

suboptimal  estimators  for  the  S system  using 
similar  techniques. 

2 

In  our  discussion  of  estimation  on  S , %re 


exp(-Jl(i4-l)(t-tQ)v)  (32) 

where  ~i***<)w*<i>*  s^erioal  har- 

monics of  degree  A (defined  in  (A.2),(A.3)).  The 
function  p(8,g.t;  U the  solution  to  (31) 

with  Initial  eondltlon  equal  to  the  singular  distri- 
bution concentrated  at  (8, q|  - (8Q.eg).  Also. 

Grenander  [20]  defines  a Gaussian  (normal)  density 
to  be  the  solution  of  (29)  for  some  t. 

The  folded  normal  density  F(8;1|.v)  used  by 
Wlllsky  as  an  assumed  density  approximation  for 
the  phase  tracking  problem  is  indeed  a nonnal  den- 
sity on  8^  in  the  sense  of  Grenander  [4].  Motiva- 
ted by  the  success  of  Willsky's  suboptimal  filter. 
\ve  will  design  suboptimal  estimators  for  bilinear 
2 

systems  on  S by  employing  a normal  assusMd 
conditional  density  of  the  form 


i«0  m*-l  * * 

• axp(-i(i>l)v(t))  (33) 

where  r(t) . X(t)  and  y{t)  ara  paramaters  of  tha  den- 
sity which  are  to  ba  estimated.  In  other  words, 
the  "generalized  Fourier  coefficients" 

2n  n 

c^^(t)  «y  J y*^(e(t).9(t))p(«.o,t)sin  tdld(p 
- E*ty^J»(0,»(t))l  (34) 


(N4-l-m)i  2N4-3 
where  a ■ — — . 

(Ni-14-m)!  4a 

Tha  decomposition  of  (A. 5)  of  hooKsganaous 

polynomials  of  degree  n (restricted  to  8^)  in  terms 
of  the  spherical  harmonics  of  degree  a n impllas 

the  existence  of  a nonsingular  matrix  P such  that 

II 

* ’^n-2^’^ 

where  Y^(x)  is  the  (2  i-t-l) -vector  whoso  components 

are  the  spherical  harmonics  (Y^  , -icmdi)  of 

sm 


degree  l and  4 is  zero  or  one  depending  on  whether 
n is  even  or  odd.  Specializing  this  result  to  the 
present  approximation  implies  the  existence  of  a 
nonsingular  matrix  P„  , such  Aat 


4(4*l)Y(t) 


In  order  to  assume  the  existence  of  the  conditional 
density,  it  is  sufficient  to  assume  a technical 
controllability  condition  [13], 


Thus  x“^  •'(t|tj  can  be  computed  from  {e^j 

-(N+1)  am  XN+l)  and  x^"*^(t|0.  The  optimal 
estimator  (27)  is  truncated  by  sulmtituting  this 

approximation  for  x^^^^(t(t)  into  dm  equation  for 


In  order  to  truncate  the  optimal  estimator 

(27)  after  the  equation  (for  x^^(t|t))  using  the 
assumed  density  (33) , we  must  compute 

E*[x^^(t)x*(t)l,  or  equivalently,  x^'*’*^(t|t),  in 

terms  of  x^^^(t|t),  p-  1,2,....N.  However,  if 
x(t|t)  is  known,  so  are  c^g(t)  and  • 

simple  computation  yields 


We  note  that  one  can  show  that 

•(t)  Ay|i(t|t)ri  *1 


and  this  quantity  can  be  used  as  a measure  of 
our  confidence  in  our  estimate.  Speciftoaiiy,  if 
x(t|t)  satisfies  the  assumed  density  (33), 


e(t)  - |x(t|t)|!  - e’^** 

so  Y ■ 0 (zero  "variance")  Impliee  1,  and 
V > • (InUnite  "variance”)  Impitea  sr  ■ 0 (sea 

for  the  enalog) . 


If  Cj|(t)  - 0,  frten  the  density  is  Independent  of 
X(t);  otherwise. 


Example;  Suppose  that  we  truncate  the 
2 

optimal  S estimator  (27)  after  M ■ 1— i.e. , wa 

approximate  x^^^(t|t)  using  tha  abova  approxima- 
tion. The  resulting  suboptimal  estimator  is  (for 
Q(t)  - I) 

dx(t|t)  - [Ag  4-  |£A^]x(t|t)dt 


axp(2U(l))  - c*,(Vc,,(t) 


(l«t) , X(t))exp(-(N>  1)(N4-2)  v(t)) 


♦ P(t)H'(DR~Nt)]dZ2(t)  - R(t)x(t|t)dd  (45) 
v/here  tha  "covariance"  matrix  P(Q  is  given  by 


(N>1)(N>2) 


4 


PjlW  * Xj(t|t)(|||x(t|t)||  - 1) 


for  1 /i.  1 k,  )/ k.  and 

Pjj(t)  - Xj(t|t)ij(t|t)(lli(t|l)||  - 1)  (47) 

for  1 y ).  NoUca  that,  from  (44),  |]x(t|t)|I  > 1 Im- 
pUes  t)>at  the  “variance"  v(t)  » 0;  In  fact,  if 
i'x(t|t)!1  « 1,  wo  aee  from  (46)-(47)  that  the  co- 
variance  matrix  P(t)  la  identically  aero.  Thus  If 
IjxU  |t)||  > 1,  this  first  order  suboptlmal  filter  as- 
sumes that  it  has  perfect  knowledge  of  x(t)  an* 
disregards  the  measurements. 

5.  Conclusion 

In  this  paper  wo  have  considered  the  esti- 
mation problem  iat  bilinear  stochastic  systems 
evolving  on  spheres.  The  techniques  of  harmonic 
analysis  have  been  applied  to  the  design  of  sub- 
optimal  estimators  for  such  systems.  The  per- 
formance of  these  estimators  is  prssentiy  being 
investigated  by  tteaiM  of  computer  simulation,  and 


(1)  the  restrlctloa  te-S^-el^he  subspaee  of  ^ 

2 f 9 

which  is  orthogonal  to  the  subspaoe  ((x'ei^exp 



2 

(2)  the  eigenspace  of  the  Laplacian  4 2 on  S 

8 

with  eigenvalue  -1(14-1),  where  [9) 

2 

V " sin*  a**'*"*^**  sin* 

Property  (1)  implies  Oat  each  f€  8^  has  a unique 
expansion  ^ 


where  ^j|-2)^^l-2)  [tl  is  dm  largest  Integer 

* t (Brockett  [9]  also  discusses  this  point) . One 
can  show  [22,  p.  109]  that  the  span  of  [H^, 

1-1 ,2 , . . . ) is  dense  In  the  space  of  continuous 
2 2 

ftinctlons  on  8 and  in  L^(8  ) , 1 e p < -. 


the  results  will  be  presented  in  a future  paper. 


Appendix  - Harmonic  Atwlysls  on  8 

Any  point  (Xj,X2,x^)  on  the  2-sphere 
2 3 

8 ” (x*  It  |x'x  - 1 } can  be  expressed  in  the  polar 
coordinates  (*,g),  where  0e*sn,  Oeg  s2n, 
by  defining 


X,  - cos*;  X,  - sin*cos9;  x.  - sin* sing 


Harmonic  aimlysis  on  S is  studied  In  terms  of  the 
spherical  harmonics  [6], [22], [23];  the  normalised 

2 

spherical  harmonica  of  degree  1 on  8 are  defined 
by  [6] 

Yi^(*,9)  ( 1)  |(2+ni)I  4n 


P.„(cos6)e 


for  1-0,1,...  and  m - 0,1,...,1,  where  P,_ 

Ain 

(cos  •)  are  the  associated  Legendre  functions  and 
* denotes  complex  cortlugate. 


Let  ^ denote  the  apace  of  homogeneous 


polynomials  of  degree  1 on  It  (l.e. , 
ficxj.cxj.cxj)  - o^f(Xj,X2,Xj)).  Then  the  space 
H^  of  spherical  harmonics  of  degree  1 on  can 
be  characterised  In  the  following  equivalent  ways; 
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